# ERP数据线性混合效应模型分析
# ERP Data Linear Mixed Effects Model Analysis

# Clean workspace
rm(list = ls())  # Remove all objects
graphics.off()    # Close all graphics devices
cat("\014")      # Clear console

# Install required packages if not already installed
if (!require(tidyverse)) install.packages("tidyverse")
if (!require(readxl)) install.packages("readxl")
if (!require(writexl)) install.packages("writexl")
if (!require(lme4)) install.packages("lme4")
if (!require(lmerTest)) install.packages("lmerTest")
if (!require(emmeans)) install.packages("emmeans")
if (!require(ggplot2)) install.packages("ggplot2")
if (!require(DHARMa)) install.packages("DHARMa")

# Load required libraries
library(tidyverse)  # for data manipulation
library(readxl)     # for reading Excel files
library(writexl)    # for writing Excel files
library(lme4)       # for mixed effects models
library(lmerTest)   # for p-values in mixed models
library(emmeans)    # for estimated marginal means
library(ggplot2)    # for visualization
library(DHARMa)     # for residual diagnostics

# Function to load ERP data
load_erp_data <- function(file_path, sheet_number = 1) {
  # Load data
  data <- read_excel(file_path, sheet = sheet_number)
  
  # Convert to tibble and clean column names
  data <- as_tibble(data)
  names(data) <- gsub(" ", "_", names(data))
  
  # Ensure proper column types
  data <- data %>%
    mutate(
      Subject = as.factor(Subject),
      Semantic_Relatedness = as.factor(Semantic_Relatedness),
      Priming_Word_Type = as.factor(Priming_Word_Type),
      Group = as.factor(Group),
      Channel = as.factor(Channel)
    )
  
  return(data)
}

# Set input file path
input_file <- "source_data/Ch_Ja_Lex_ERP_2Group_for R.xlsx"

# Create result folder name from input file name
result_folder <- tools::file_path_sans_ext(basename(input_file))
result_folder <- gsub("_for R", "", result_folder)
result_path <- file.path("results", result_folder, "mixed_effects_analysis")

# Create result folder
if (!dir.exists(result_path)) {
  dir.create(result_path, recursive = TRUE)
  print(paste("Created result folder:", result_path))
}

# Load ERP data
print("Loading ERP data from sheet 3...")
erp_data <- load_erp_data(input_file, sheet_number = 3)

# Define subjects to exclude
excluded_subjects <- c(75, 77, 79, 80, 81, 83, 84, 85)

# Filter out excluded subjects
print("Filtering out excluded subjects...")
erp_data <- erp_data %>%
  filter(!(Subject %in% excluded_subjects))

# Print summary of remaining subjects
print(paste("Number of subjects after filtering:", length(unique(erp_data$Subject))))
cat("Excluded subjects:", paste(excluded_subjects, collapse = ", "), "\n")

# Define electrodes for different time windows
electrodes_190_250 <- c("FP1", "FP2", "Fz", "F3", "F4", "FCz", "FC3", "FC4")
electrodes_other <- c("Cz", "C3", "C4", "FCz", "FC1", "FC2", "CPz", "CP1", "CP2")

# Function to analyze one ERP time window
analyze_erp_window <- function(data, erp_var, result_path) {
  print(paste("\nAnalyzing", erp_var, "..."))
  
  # Select appropriate electrodes based on time window
  if (erp_var == "ERP_190_250") {
    print("Using frontal electrodes for 190-250ms time window...")
    analysis_data <- data %>%
      filter(Channel %in% electrodes_190_250)
    used_electrodes <- electrodes_190_250
  } else {
    print("Using central electrodes for other time windows...")
    analysis_data <- data %>%
      filter(Channel %in% electrodes_other)
    used_electrodes <- electrodes_other
  }
  
  # Print electrode information
  print(paste("Number of electrodes used:", length(used_electrodes)))
  print("Selected electrodes:")
  print(used_electrodes)
  
  # Set up model formula
  model_formula <- as.formula(paste0(
    erp_var, " ~ (Semantic_Relatedness + Priming_Word_Type + Group)^2 + ",
    "Semantic_Relatedness:Priming_Word_Type:Group + ",
    "(1 + Semantic_Relatedness + Priming_Word_Type | Subject) + ",
    "(1 + Semantic_Relatedness | Group:Channel)"
  ))
  
  # Fit the model
  print("Fitting linear mixed-effects model...")
  mixed_model <- lmer(model_formula, data = analysis_data)
  print("Model fitting completed")
  
  # Print model summary
  print("\n============= Model Summary =============")
  model_summary <- summary(mixed_model)
  print(model_summary)
  
  # Extract fixed effects and their significance
  fixed_effects <- fixef(mixed_model)
  fixed_effects_se <- sqrt(diag(vcov(mixed_model)))
  fixed_effects_t <- fixed_effects / fixed_effects_se
  fixed_effects_p <- 2 * (1 - pnorm(abs(fixed_effects_t)))
  
  # Add significance markers function
  get_significance_markers <- function(p_value) {
    if (p_value < 0.001) return(" ***")
    else if (p_value < 0.01) return(" **")
    else if (p_value < 0.05) return(" *")
    else if (p_value < 0.1) return(" .")
    else return("")
  }
  
  # Create fixed effects table with significance markers
  fixed_effects_table <- data.frame(
    Estimate = fixed_effects,
    SE = fixed_effects_se,
    t_value = fixed_effects_t,
    p_value = paste0(round(fixed_effects_p, 4), sapply(fixed_effects_p, get_significance_markers))  # Add markers directly to p-values
  )
  
  # Create results report file
  report_file <- file.path(result_path, "ERP_分析报告.txt")
  
  # Write report header
  cat(paste0(
    "\n\n", erp_var, " 时间窗口分析报告\n",
    "================================\n\n",
    "1. 模型概述\n",
    "   时间窗口：", erp_var, "\n",
    "   观测数：", nrow(analysis_data), "\n",
    "   被试数：", length(unique(analysis_data$Subject)), "\n",
    "   通道数：", length(unique(analysis_data$Channel)), "\n",
    "   使用的电极点：", paste(used_electrodes, collapse = ", "), "\n\n"
  ), file = report_file, append = TRUE, fileEncoding = "UTF-8")
  
  # Write fixed effects results
  cat("\n2. 固定效应结果\n", file = report_file, append = TRUE)
  capture.output(
    {
      print(fixed_effects_table, digits = 4)
      cat("\nSignificance codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n")
    },
    file = report_file,
    append = TRUE
  )
  
  # Add random effects results
  cat("\n3. 随机效应结果\n", file = report_file, append = TRUE)
  random_effects <- VarCorr(mixed_model)
  capture.output(
    print(random_effects, comp = c("Variance", "Std.Dev.")),
    file = report_file,
    append = TRUE
  )
  
  # Perform simple effects analysis
  print("\nPerforming simple effects analysis...")
  
  # Simple effects by Group
  emm_group <- emmeans(mixed_model, pairwise ~ Priming_Word_Type|Group)
  # Simple effects by Semantic Relatedness
  emm_relatedness <- emmeans(mixed_model, pairwise ~ Group|Semantic_Relatedness)
  # Simple effects by Priming Word Type
  emm_priming <- emmeans(mixed_model, pairwise ~ Group|Priming_Word_Type)
  
  # Add detailed analysis of Semantic_Relatedness × Priming_Word_Type interaction for each group
  print("\nAnalyzing Semantic_Relatedness × Priming_Word_Type interaction by group...")
  emm_group_interaction <- emmeans(mixed_model, pairwise ~ Semantic_Relatedness:Priming_Word_Type|Group)
  
  # Add simple effects to report
  cat("\n4. 简单效应分析\n", file = report_file, append = TRUE)
  capture.output(
    {
      cat("\n4.1 按组别的简单效应：\n")
      print(emm_group)
      cat("\n4.2 按语义关联性的简单效应：\n")
      print(emm_relatedness)
      cat("\n4.3 按启动词类型的简单效应：\n")
      print(emm_priming)
      cat("\n4.4 各组内语义关联性与启动词类型的交互效应：\n")
      print(emm_group_interaction)
    },
    file = report_file,
    append = TRUE
  )
  
  # Add interpretation with significance markers
  cat(paste0(
    "\n5. 结果解释\n",
    "   5.1 主效应分析\n",
    "   - Group主效应：", ifelse(fixed_effects_p["Group1"] < 0.05, 
                             paste0("显著 (p = ", round(fixed_effects_p["Group1"], 4), get_significance_markers(fixed_effects_p["Group1"])),
                             paste0("不显著 (p = ", round(fixed_effects_p["Group1"], 4), ")")), "\n",
    "   - Semantic_Relatedness主效应：", ifelse(fixed_effects_p["Semantic_Relatedness1"] < 0.05,
                                           paste0("显著 (p = ", round(fixed_effects_p["Semantic_Relatedness1"], 4), get_significance_markers(fixed_effects_p["Semantic_Relatedness1"])),
                                           paste0("不显著 (p = ", round(fixed_effects_p["Semantic_Relatedness1"], 4), ")")), "\n",
    "   - Priming_Word_Type主效应：", ifelse(fixed_effects_p["Priming_Word_Type1"] < 0.05,
                                        paste0("显著 (p = ", round(fixed_effects_p["Priming_Word_Type1"], 4), get_significance_markers(fixed_effects_p["Priming_Word_Type1"])),
                                        paste0("不显著 (p = ", round(fixed_effects_p["Priming_Word_Type1"], 4), ")")), "\n\n",
    
    "   5.2 交互作用分析\n",
    "   - Group × Semantic_Relatedness：", ifelse(fixed_effects_p["Group1:Semantic_Relatedness1"] < 0.05,
                                             paste0("显著 (p = ", round(fixed_effects_p["Group1:Semantic_Relatedness1"], 4), get_significance_markers(fixed_effects_p["Group1:Semantic_Relatedness1"])),
                                             paste0("不显著 (p = ", round(fixed_effects_p["Group1:Semantic_Relatedness1"], 4), ")")), "\n",
    "   - Group × Priming_Word_Type：", ifelse(fixed_effects_p["Group1:Priming_Word_Type1"] < 0.05,
                                          paste0("显著 (p = ", round(fixed_effects_p["Group1:Priming_Word_Type1"], 4), get_significance_markers(fixed_effects_p["Group1:Priming_Word_Type1"])),
                                          paste0("不显著 (p = ", round(fixed_effects_p["Group1:Priming_Word_Type1"], 4), ")")), "\n",
    "   - Semantic_Relatedness × Priming_Word_Type：", ifelse(fixed_effects_p["Semantic_Relatedness1:Priming_Word_Type1"] < 0.05,
                                                        paste0("显著 (p = ", round(fixed_effects_p["Semantic_Relatedness1:Priming_Word_Type1"], 4), get_significance_markers(fixed_effects_p["Semantic_Relatedness1:Priming_Word_Type1"])),
                                                        paste0("不显著 (p = ", round(fixed_effects_p["Semantic_Relatedness1:Priming_Word_Type1"], 4), ")")), "\n",
    "   - 三阶交互：", ifelse(fixed_effects_p["Group1:Semantic_Relatedness1:Priming_Word_Type1"] < 0.05,
                          paste0("显著 (p = ", round(fixed_effects_p["Group1:Semantic_Relatedness1:Priming_Word_Type1"], 4), get_significance_markers(fixed_effects_p["Group1:Semantic_Relatedness1:Priming_Word_Type1"])),
                          paste0("不显著 (p = ", round(fixed_effects_p["Group1:Semantic_Relatedness1:Priming_Word_Type1"], 4), ")")), "\n\n",
    
    "   5.3 各组内语义关联性与启动词类型的交互效应\n",
    "   详细结果请参见4.4节\n\n",
    
    "6. 随机效应解释\n",
    "   - Subject随机效应：表示个体间ERP波幅的变异程度\n",
    "   - Channel随机效应：表示不同电极通道间的变异程度\n\n",
    
    "7. 总结\n",
    "   根据以上分析结果，主要发现：\n",
    "   - ", ifelse(any(fixed_effects_p[c("Group1", "Semantic_Relatedness1", "Priming_Word_Type1")] < 0.05),
                 "存在显著的主效应：\n     ",
                 "未发现显著的主效应\n"), 
    ifelse(fixed_effects_p["Group1"] < 0.05, "Group效应显著\n     ", ""),
    ifelse(fixed_effects_p["Semantic_Relatedness1"] < 0.05, "语义关联性效应显著\n     ", ""),
    ifelse(fixed_effects_p["Priming_Word_Type1"] < 0.05, "启动词类型效应显著\n     ", ""),
    "   - ", ifelse(any(fixed_effects_p[grep(":", names(fixed_effects_p))] < 0.05),
                 "存在显著的交互作用：\n     ",
                 "未发现显著的交互作用\n"),
    "\n注：所有效应的具体方向和大小请参考上述详细统计结果。\n\n"
  ), file = report_file, append = TRUE, fileEncoding = "UTF-8")
  
  # Return model and summary for potential further analysis
  return(list(
    model = mixed_model,
    summary = model_summary,
    fixed_effects = fixed_effects_table
  ))
}

# Analyze each ERP time window
results_list <- list()
for (window in c("ERP_190_250", "ERP_330_480", "ERP_500_700")) {
  results_list[[window]] <- analyze_erp_window(erp_data, window, result_path)
}

# Add final timestamp to the report
cat(paste0("\n\n分析日期：", format(Sys.Date(), "%Y年%m月%d日"), "\n"),
    file = file.path(result_path, "ERP_分析报告.txt"),
    append = TRUE, fileEncoding = "UTF-8")

print("\nAnalysis completed for all ERP time windows.")
print(paste("Results saved in:", result_path)) 